#***************************************************************#
#                                                               #
#                 Replication file (PART 1) for:                #
#                                                               #
#         "The Effect of Anger Appeals on the Support for       #
#                     Secessionist Parties"                     #
#                        Barceló, Joan                          #
#                  Journal of Politics (JOP)                    #
#               Last updated: 10 December 2022                  #
#                                                               #
#****************************************************************

#################
#               #
# Load packages #
#               #
#################

library(rstudioapi)
library(foreign)
library(dplyr)
library(tidyverse)
library(htmlTable)
library(compareDF)
library(AER)
library(ivreg)
library(reshape2)
library(stargazer)
library(outliers)

#############
# Load Data #
#############

current_path <- getActiveDocumentContext()$path 
setwd(dirname(current_path))

FieldExperimentData <- readRDS('FieldExperimentData.RDS')

###########
# Table 1 #
###########

#Vote for pro-independence parties in 2015
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$proindependence2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$proindependence2015, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$proindependence2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$proindependence2015, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$proindependence2015)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$proindependence2015)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$proindependence2015)
orto_independence15 <- aov(proindependence2015 ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_independence15)

#Vote for ambivalent parties in 2015
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$ambivalent2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$ambivalent2015, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$ambivalent2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$ambivalent2015, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$ambivalent2015)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$ambivalent2015)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$ambivalent2015)
orto_ambivalent15 <- aov(ambivalent2015 ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_ambivalent15)

#Turnout in 2015
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$turnout2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$Turnout2015, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$turnout2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$Turnout2015, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$turnout2015)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$turnout2015)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$turnout2015)
orto_turnout15 <- aov(turnout2015 ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_turnout15)

#Number of eligible voters in 2015
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$census2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$census2015, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$census2015,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$census2015, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$census2015)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$census2015)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$census2015)
orto_census15 <- aov(census2015 ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_census15)

#Distance to province capital
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$distance_prov_capital,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$distance_prov_capital, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$distance_prov_capital,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$distance_prov_capital, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$distance_prov_capital)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$distance_prov_capital)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$distance_prov_capital)
orto_capdist <- aov(distance_prov_capital ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_capdist)

#Population size
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$population_size,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$population_size, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$population_size,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$population_size, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$population_size)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$population_size)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$population_size)
orto_popsize <- aov(population_size ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_popsize)

#Population density
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$population_density,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$population_density, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$population_density,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$population_density, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$population_density)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$population_density)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$population_density)
orto_density<- aov(population_density ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_density)

#Share of Catalan-born population
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$catalanborn_share,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$catalanborn_share, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$catalanborn_share,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$catalanborn_share, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$catalanborn_share)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$catalanborn_share)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$catalanborn_share)
orto_catalanborn <- aov(catalanborn_share ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_catalanborn)

#Share with prior experience with police brutality
t.test(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$police_violence,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$police_violence, alternative = "two.sided", var.equal = FALSE)
t.test(FieldExperimentData[which(FieldExperimentData$control == 1),]$police_violence,FieldExperimentData[which(FieldExperimentData$placebo == 1),]$police_violence, alternative = "two.sided", var.equal = FALSE)
sd(FieldExperimentData[which(FieldExperimentData$treatment == 1),]$police_violence)
sd(FieldExperimentData[which(FieldExperimentData$placebo == 1),]$police_violence)
sd(FieldExperimentData[which(FieldExperimentData$control == 1),]$police_violence)
orto_violence <- aov(police_violence ~ letter_condition, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_violence)

#Orthogonality tests

orto_full_placebo <- lm(placebo ~ proindependence2015 + ambivalent2015 + turnout2015 + police_violence + census2015 + distance_prov_capital + population_density + population_size + catalanborn_share, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '2'),])
summary(orto_full_placebo)

orto_full_treatment <- lm(treatment ~ proindependence2015 + ambivalent2015 + turnout2015 + police_violence + census2015 + distance_prov_capital + population_density + population_size + catalanborn_share, data = FieldExperimentData[which(FieldExperimentData$letter_condition == '1' |  FieldExperimentData$letter_condition == '3'),])
summary(orto_full_treatment)

##################################
# Figure 4 and online Appendix G #
##################################

#Plot Main Effects: ITT

summary(mod_indep0 <- lm(proindependence2017 ~ as.factor(placebo) + as.factor(treatment) + proindependence2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
res_indep0 <- summary(mod_indep0)
pt(coef(res_indep0)[, 3], mod_indep0$df, lower = FALSE)

summary(mod_indep1 <- lm(proindependence2017 ~ as.factor(treatment) + proindependence2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
res_indep1 <- summary(mod_indep1)
pt(coef(res_indep1)[, 3], mod_indep1$df, lower = FALSE)

#Plot Main Effects: CACE
summary(cace_indep_wp <- ivreg(proindependence2017 ~ treat_rate + placebo_rate + proindependence2015 | as.factor(treatment) + as.factor(placebo) + proindependence2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
summary(cace_indep <- ivreg(proindependence2017 ~ treat_rate + proindependence2015 | as.factor(treatment) + proindependence2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))

# Estimates into temporary data.frames:
model1Frame <- data.frame(Variable = c(' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter'),
                          Coefficient = c(0, summary(mod_indep0)$coefficients[2,1], summary(mod_indep0)$coefficients[3,1], 0, 0, summary(mod_indep1)$coefficients[2,1], 0, summary(cace_indep_wp)$coefficients[3,1], summary(cace_indep_wp)$coefficients[2,1], 0, 0, summary(cace_indep)$coefficients[2,1]),
                          SE = c(0, summary(mod_indep0)$coefficients[2,2], summary(mod_indep0)$coefficients[3,2], 0, 0, summary(mod_indep1)$coefficients[2,2], 0, summary(cace_indep_wp)$coefficients[3,2], summary(cace_indep_wp)$coefficients[2,2], 0, 0, summary(cace_indep)$coefficients[2,2]),
                          Model = c('1. ITT (Treatment/Placebo vs. Control)', '1. ITT (Treatment/Placebo vs. Control)', '1. ITT (Treatment/Placebo vs. Control)', 
                                    '2. ITT (Treatment vs. Any Other Condition)', '2. ITT (Treatment vs. Any Other Condition)', '2. ITT (Treatment vs. Any Other Condition)', 
                                    ' 3. CACE (Treatment/Placebo vs. Control)', ' 3. CACE (Treatment/Placebo vs. Control)', ' 3. CACE (Treatment/Placebo vs. Control)', 
                                    ' 4. CACE (Treatment vs. Any Other Condition)', ' 4. CACE (Treatment vs. Any Other Condition)', ' 4. CACE (Treatment vs. Any Other Condition)'),
                          show0 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show1 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show2 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show3 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE))

allModelFrame <- data.frame(rbind(model1Frame))  # etc.

# Width of CI
interval1 <- -qnorm((1-0.9))  # 90% multiplier, one-tailed
interval2 <- -qnorm((1-0.95))  # 95% multiplier, one-tailed

# Plot

allModelFrame$Variable <- factor(allModelFrame$Variable)
allModelFrame$Model <- factor(allModelFrame$Model)
allModelFrame$Model <- factor(allModelFrame$Model,levels(allModelFrame$Model)[c(3,4,1,2)])

zp1 <- ggplot(allModelFrame, aes(group=Model, colour=Model, alpha = show0))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + labs(x="",y="Estimated Effects on \n Pro-independence Parties Vote Shares")
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1.5, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2, ymax = Coefficient + SE*interval2),
                             lwd = 0.75, position = position_dodge(width = 1/2),
                             shape = 19)
zp1 <- zp1 + theme_bw() + theme(axis.text=element_text(size=14),
                                axis.title=element_text(size=18)) 
zp1 <- zp1 + ggtitle("") + scale_alpha_identity()
print(zp1)


###################################
# Figure 5A and online Appendix I #
###################################

summary(mod_ambivalent0 <- lm(ambivalent2017 ~ as.factor(placebo) + as.factor(treatment) + ambivalent2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
t_stats0 = (mod_ambivalent0$coefficients - 0)/coef(summary(mod_ambivalent0))[, 2]
p_values0 = pt(t_stats0, mod_ambivalent0$df.residual)

summary(mod_ambivalent1 <- lm(ambivalent2017 ~ as.factor(treatment) + ambivalent2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
t_stats1 = (mod_ambivalent1$coefficients - 0)/coef(summary(mod_ambivalent1))[, 2]
p_values1 = pt(t_stats1, mod_ambivalent1$df.residual)

cace_ambivalent_wp <- ivreg(ambivalent2017 ~ treat_rate + placebo_rate + ambivalent2015 | as.factor(treatment) + as.factor(placebo) + ambivalent2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,])

cace_ambivalent <- ivreg(ambivalent2017 ~ treat_rate + ambivalent2015 | as.factor(treatment) + ambivalent2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,])

modelNFrame <- data.frame(Variable = c(' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter'),
                          Coefficient = c(0, summary(mod_ambivalent0)$coefficients[2,1], summary(mod_ambivalent0)$coefficients[3,1], 0, 0, summary(mod_ambivalent1)$coefficients[2,1], 0, summary(cace_ambivalent_wp)$coefficients[3,1], summary(cace_ambivalent_wp)$coefficients[2,1], 0, 0, summary(cace_ambivalent)$coefficients[2,1]),
                          SE = c(0, summary(mod_ambivalent0)$coefficients[2,2], summary(mod_ambivalent0)$coefficients[3,2], 0, 0, summary(mod_ambivalent1)$coefficients[2,2], 0, summary(cace_ambivalent_wp)$coefficients[3,2], summary(cace_ambivalent_wp)$coefficients[2,2], 0, 0, summary(cace_ambivalent)$coefficients[2,2]),
                          Model = c('1. ITT (Treatment/Placebo vs. Control)', '1. ITT (Treatment/Placebo vs. Control)', '1. ITT (Treatment/Placebo vs. Control)', 
                                    '2. ITT (Treatment vs. Any Other Condition)', '2. ITT (Treatment vs. Any Other Condition)', '2. ITT (Treatment vs. Any Other Condition)', 
                                    ' 3. CACE (Treatment/Placebo vs. Control)', ' 3. CACE (Treatment/Placebo vs. Control)', ' 3. CACE (Treatment/Placebo vs. Control)', 
                                    ' 4. CACE (Treatment vs. Any Other Condition)', ' 4. CACE (Treatment vs. Any Other Condition)', ' 4. CACE (Treatment vs. Any Other Condition)'),
                          show0 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show1 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show2 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show3 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE))

allModelFrameN <- data.frame(rbind(modelNFrame))  # etc.

# Width of CI
interval1 <- -qnorm((1-0.9))  # 90% multiplier, one-tailed
interval2 <- -qnorm((1-0.95))  # 95% multiplier, one-tailed

# Plot

allModelFrameN$Variable <- factor(allModelFrameN$Variable)
allModelFrameN$Model <- factor(allModelFrameN$Model)
allModelFrameN$Model <- factor(allModelFrameN$Model,levels(allModelFrameN$Model)[c(3,4,1,2)])

np1 <- ggplot(allModelFrameN, aes(group=Model, colour=Model, alpha = show0))
np1 <- np1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + labs(x="",y="Estimated Effects on \n Voting for Ambivalent parties")
np1 <- np1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1.5, position = position_dodge(width = 1/2))
np1 <- np1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2, ymax = Coefficient + SE*interval2),
                             lwd = 0.75, position = position_dodge(width = 1/2),
                             shape = 19)
np1 <- np1 + theme_bw() + theme(axis.text=element_text(size=14),
                                axis.title=element_text(size=18)) 
np1 <- np1 + ggtitle("") + scale_alpha_identity()
print(np1)


###################################
# Figure 5B and online Appendix I #
###################################

#Plot Effects on Unionist Parties: ITT
summary(mod_union0 <- lm(prounion2017 ~ as.factor(placebo) + as.factor(treatment) + prounion2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
t_stats0 = (mod_union0$coefficients - 0)/coef(summary(mod_union0))[, 2]
p_values0 = 1*pt(t_stats0, mod_union0$df.residual, lower = FALSE)

summary(mod_union1 <- lm(prounion2017 ~ as.factor(treatment) + prounion2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
t_stats1 = (mod_union1$coefficients - 0)/coef(summary(mod_union1))[, 2]
p_values1 = 1*pt(t_stats1, mod_union1$df.residual, lower = FALSE)

summary(cace_union <- ivreg(prounion2017 ~ treat_rate + prounion2015 | as.factor(treatment) + prounion2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))

summary(cace_union_wp <- ivreg(prounion2017 ~ treat_rate + placebo_rate + prounion2015 | as.factor(treatment) + as.factor(placebo) + prounion2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))

modelUFrame <- data.frame(Variable = c(' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter', ' No letter', 'Emotionally\n Neutral \nletter', 'Violence/Anger \nletter'),
                          Coefficient = c(0, summary(mod_union0)$coefficients[2,1], summary(mod_union0)$coefficients[3,1], 0, 0, summary(mod_union1)$coefficients[2,1], 0, summary(cace_union_wp)$coefficients[3,1], summary(cace_union_wp)$coefficients[2,1], 0, 0, summary(cace_union)$coefficients[2,1]),
                          SE = c(0, summary(mod_union0)$coefficients[2,2], summary(mod_union0)$coefficients[3,2], 0, 0, summary(mod_union1)$coefficients[2,2], 0, summary(cace_union_wp)$coefficients[3,2], summary(cace_union_wp)$coefficients[2,2], 0, 0, summary(cace_union)$coefficients[2,2]),
                          Model = c('1. ITT (Treatment/Placebo vs. Control)', '1. ITT (Treatment/Placebo vs. Control)', '1. ITT (Treatment/Placebo vs. Control)', 
                                    '2. ITT (Treatment vs. Any Other Condition)', '2. ITT (Treatment vs. Any Other Condition)', '2. ITT (Treatment vs. Any Other Condition)', 
                                    ' 3. CACE (Treatment/Placebo vs. Control)', ' 3. CACE (Treatment/Placebo vs. Control)', ' 3. CACE (Treatment/Placebo vs. Control)', 
                                    ' 4. CACE (Treatment vs. Any Other Condition)', ' 4. CACE (Treatment vs. Any Other Condition)', ' 4. CACE (Treatment vs. Any Other Condition)'),
                          show0 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show1 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show2 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE),
                          show3 = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE))

allModelFrameU <- data.frame(rbind(modelUFrame))  # etc.

# Width of CI
interval1 <- -qnorm((1-0.9))  # 90% multiplier, one-tailed
interval2 <- -qnorm((1-0.95))  # 95% multiplier, one-tailed

# Plot
allModelFrameU$Variable <- factor(allModelFrameU$Variable)
allModelFrameU$Model <- factor(allModelFrameU$Model)
allModelFrameU$Model <- factor(allModelFrameU$Model,levels(allModelFrameU$Model)[c(3,4,1,2)])

up1 <- ggplot(allModelFrameU, aes(group=Model, colour=Model, alpha = show0))
up1 <- up1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + labs(x="",y="Estimated Effects on \n Voting for Unionist parties")
up1 <- up1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1.5, position = position_dodge(width = 1/2))
up1 <- up1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2, ymax = Coefficient + SE*interval2),
                             lwd = 0.75, position = position_dodge(width = 1/2),
                             shape = 19)
up1 <- up1 + theme_bw() + theme(axis.text=element_text(size=14),
                                axis.title=element_text(size=18)) 
up1 <- up1 + ggtitle("") + scale_alpha_identity()
print(up1)

######################
# Online Appendix F  #
######################

summary(mode_not_reached_outcomes_out1 <- lm(not_reached ~ proindependence2015, data = FieldExperimentData[!(FieldExperimentData$letter_condition == 0 | FieldExperimentData$letter_condition == 1),]))
linearHypothesis(mode_not_reached_outcomes_out1, c("proindependence2015=0"))

summary(mode_not_reached_outcomes_out2 <- lm(not_reached ~ proindependence2015 + ambivalent2015, data = FieldExperimentData[!(FieldExperimentData$letter_condition == 0 | FieldExperimentData$letter_condition == 1),]))
linearHypothesis(mode_not_reached_outcomes_out2, c("proindependence2015=0", "ambivalent2015=0"))
linearHypothesis(mode_not_reached_outcomes_out2, c("ambivalent2015=0"))

######################
# Online Appendix H  #
######################

set.seed(123)
n = nrow(FieldExperimentData[!FieldExperimentData$letter_condition == 0,])
n_treat = nrow(FieldExperimentData[FieldExperimentData$letter_condition == 3,])
Y = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$proindependence2017 - FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$proindependence2015
compl_rate = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$compliance_rate

#ProIndependence

mean(replicate(100000, {
  # randomly assign treatment or control
  treat = rbinom(n, 1, 0.5)
  # get the observed outcome
  # calculate the difference in means
  diff_means = mean(Y[treat == 1]) - mean(Y[treat == 0])
  diff_means
}))

# Draw a randomization used for the experiment
treat_obs = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$treatment

# Randomization distribution

ri_dist_itt = replicate(100000, {
  treat_new = as.numeric(1:n %in% sample(1:n, n_treat, FALSE))
  diff_means = mean(Y[treat_new == 1]) - mean(Y[treat_new == 0])
})

diff_means_obs = round(mod_indep1$coefficients[2], 3)

p_val = round(sum(ri_dist_itt > diff_means_obs) / length(ri_dist_itt), 3)
p_val

g <- ggplot(data = data.frame(i = 1:length(ri_dist_itt), ri_dist = ri_dist_itt)) + geom_histogram(aes(ri_dist_itt, fill=I(rgb(0/256, 0/256, 0/256))), bins = 50) + xlab("ITT under Sharp Null (100,000 simulations)") + ylab("Count") + #ggtitle("Randomization Distribution of Difference in Means Under a Strict Null") 
  geom_vline(xintercept = diff_means_obs, linetype = "dashed", col = 'red') + annotate("text", x = diff_means_obs, y = 4500, label = paste0("Observed ITT: ", round(diff_means_obs, 2), "\n P-value: ", p_val), hjust = -0.1, size = 6) + theme_bw()

g + theme(axis.text=element_text(size=20), axis.title=element_text(size=22, face="bold"))

##Sharp Null CACE

Y = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$proindependence2017
Y_1 = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$proindependence2015

diff_means = matrix()
mean(replicate(10000, {
  # randomly assign treatment or control
  treat = rbinom(n, 1, 0.5)
  #Calculate compliance rate
  Z = treat*compl_rate
  diff_means = coefficients(ivreg(Y ~ Z + Y_1 | treat + Y_1, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))[2]
  diff_means
}))

#CACE: Pro-Independence
set.seed(123)

ri_dist_cace = replicate(100000, {
  treat_new = as.numeric(1:n %in% sample(1:n, n_treat, FALSE))
  Z = treat_new*compl_rate
  diff_means = coefficients(ivreg(Y ~ Z + Y_1 | treat_new + Y_1))[2]
  diff_means})

diff_means_obs = round(coefficients(ivreg(proindependence2017 ~ treat_rate + proindependence2015 | as.factor(treatment) + proindependence2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))[2], 3)

p_val = round(sum(ri_dist_cace > diff_means_obs) / length(ri_dist_cace), 3)
p_val

g <- ggplot(data = data.frame(i = 1:length(ri_dist_cace), ri_dist = ri_dist_cace)) + geom_histogram(aes(ri_dist_cace, fill=I(rgb(0/256, 0/256, 0/256))), bins = 50) + xlab("CACE under Sharp Null (100,000 simulations)") + ylab("Count") + #ggtitle("Randomization Distribution of Difference in Means Under a Strict Null") 
  geom_vline(xintercept = diff_means_obs, linetype = "dashed", col = 'red') + annotate("text", x = diff_means_obs, y = 4500, label = paste0("Observed CACE: ", round(diff_means_obs, 2), "\n P-value: ", p_val), hjust = -0.1, size = 6) + theme_bw()

g + theme(axis.text=element_text(size=20), axis.title=element_text(size=22, face="bold"))


######################
# Online Appendix J  #
######################

Y_ambivalent = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$ambivalent2017 - FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$ambivalent2015
Y_union = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$prounion2017 - FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$prounion2015

#ITT: Ambivalent
set.seed(123)
mean(replicate(100000, {
  # randomly assign treatment or control
  treat = rbinom(n, 1, 0.5)
  # get the observed outcome
  # calculate the difference in means
  diff_means = mean(Y_ambivalent[treat == 1]) - mean(Y_ambivalent[treat == 0])
  diff_means
}))

set.seed(123)
ri_dist_itt_ambivalent = replicate(100000, {
  treat_new = as.numeric(1:n %in% sample(1:n, n_treat, FALSE))
  diff_means = mean(Y_ambivalent[treat_new == 1]) - mean(Y_ambivalent[treat_new == 0])
})

diff_means_obs = round(mod_ambivalent1$coefficients[2], 3)

p_val = round(sum(ri_dist_itt_ambivalent < diff_means_obs) / length(ri_dist_itt_ambivalent), 3)
p_val

g <- ggplot(data = data.frame(i = 1:length(ri_dist_itt_ambivalent), ri_dist = ri_dist_itt_ambivalent)) + geom_histogram(aes(ri_dist_itt_ambivalent, fill=I(rgb(0/256, 0/256, 0/256))), bins = 50) + xlab("ITT under Sharp Null (100,000 simulations)") + ylab("Count") + #ggtitle("Randomization Distribution of Difference in Means Under a Strict Null") 
  geom_vline(xintercept = diff_means_obs, linetype = "dashed", col = 'red') + annotate("text", x = diff_means_obs-2, y = 4500, label = paste0("Observed ITT: ", round(diff_means_obs, 2), "\n P-value: ", p_val), hjust = -0.1, size = 6) + theme_bw()

g + theme(axis.text=element_text(size=20), axis.title=element_text(size=22, face="bold"))

#ITT: Unionist

mean(replicate(100000, {
  # randomly assign treatment or control
  treat = rbinom(n, 1, 0.5)
  # get the observed outcome
  # calculate the difference in means
  diff_means = mean(Y_union[treat == 1]) - mean(Y_union[treat == 0])
  diff_means
}))

set.seed(123)
ri_dist_itt_union = replicate(100000, {
  treat_new = as.numeric(1:n %in% sample(1:n, n_treat, FALSE))
  diff_means = mean(Y_union[treat_new == 1]) - mean(Y_union[treat_new == 0])
})

diff_means_obs = round(mod_union1$coefficients[2], 3)

p_val = round(sum(ri_dist_itt_union < diff_means_obs) / length(ri_dist_itt_union), 3)
p_val

g <- ggplot(data = data.frame(i = 1:length(ri_dist_itt_union), ri_dist = ri_dist_itt_union)) + geom_histogram(aes(ri_dist_itt_union, fill=I(rgb(0/256, 0/256, 0/256))), bins = 50) + xlab("ITT under Sharp Null (100,000 simulations)") + ylab("Count") + #ggtitle("Randomization Distribution of Difference in Means Under a Strict Null") 
  geom_vline(xintercept = diff_means_obs, linetype = "dashed", col = 'red') + annotate("text", x = diff_means_obs+0.35, y = 4500, label = paste0("Observed ITT: ", round(diff_means_obs, 2), "\n P-value: ", p_val), hjust = -0.1, size = 6) + theme_bw()

g + theme(axis.text=element_text(size=20), axis.title=element_text(size=22, face="bold"))


#CACE: Ambivalent Parties

Y_ambivalent = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$ambivalent2017
Y_ambivalent_1 = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$ambivalent2015

ri_dist_cace_ambivalent = replicate(100000, {
  treat_new = as.numeric(1:n %in% sample(1:n, n_treat, FALSE))
  Z = treat_new*compl_rate
  diff_means = coefficients(ivreg(Y_ambivalent ~ Z + Y_ambivalent_1 | treat_new + Y_ambivalent_1))[2]
  diff_means})

diff_means_obs = round(coefficients(ivreg(ambivalent2017 ~ treat_rate + ambivalent2015 | as.factor(treatment) + ambivalent2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))[2], 3)

p_val = round(sum(ri_dist_cace_ambivalent < diff_means_obs) / length(ri_dist_cace_ambivalent), 3)
p_val

g <- ggplot(data = data.frame(i = 1:length(ri_dist_cace_ambivalent), ri_dist = ri_dist_cace_ambivalent)) + geom_histogram(aes(ri_dist_cace_ambivalent, fill=I(rgb(0/256, 0/256, 0/256))), bins = 50) + xlab("CACE under Sharp Null (100,000 simulations)") + ylab("Count") + #ggtitle("Randomization Distribution of Difference in Means Under a Strict Null") 
  geom_vline(xintercept = diff_means_obs, linetype = "dashed", col = 'red') + annotate("text", x = diff_means_obs-1.7, y = 4500, label = paste0("Observed CACE: ", round(diff_means_obs, 2), "\n P-value: ", p_val), hjust = -0.1, size = 6) + theme_bw()

g + theme(axis.text=element_text(size=20), axis.title=element_text(size=22, face="bold"))

#CACE: Unionist Parties

Y_union = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$prounion2017
Y_union_1 = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]$prounion2015

ri_dist_cace_union = replicate(100000, {
  treat_new = as.numeric(1:n %in% sample(1:n, n_treat, FALSE))
  Z = treat_new*compl_rate
  diff_means = coefficients(ivreg(Y_union ~ Z + Y_union_1 | treat_new + Y_union_1))[2]
  diff_means})

diff_means_obs = round(coefficients(cace_union)[2], 3)

p_val = round(sum(ri_dist_cace_union < diff_means_obs) / length(ri_dist_cace_union), 3)
p_val

g <- ggplot(data = data.frame(i = 1:length(ri_dist_cace_union), ri_dist = ri_dist_cace_union)) + geom_histogram(aes(ri_dist_cace_union, fill=I(rgb(0/256, 0/256, 0/256))), bins = 50) + xlab("CACE under Sharp Null (100,000 simulations)") + ylab("Count") + #ggtitle("Randomization Distribution of Difference in Means Under a Strict Null") 
  geom_vline(xintercept = diff_means_obs, linetype = "dashed", col = 'red') + annotate("text", x = diff_means_obs+0.6, y = 4500, label = paste0("Observed CACE: ", round(diff_means_obs, 2), "\n P-value: ", p_val), hjust = -0.1, size = 6) + theme_bw()

g + theme(axis.text=element_text(size=20), axis.title=element_text(size=22, face="bold"))

#############################################
# In-text information and online Appendix K #
#############################################
#############
# Table K.1 #
#############

#Column 1
summary(mod_turnout <- lm(turnout2017 ~ as.factor(placebo) + as.factor(treatment) + turnout2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
#Column 2
summary(mod_turnout1 <- lm(turnout2017 ~ as.factor(treatment) + turnout2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,]))
#Column 3
cace_turnout <- ivreg(turnout2017 ~ treat_rate + placebo_rate + turnout2015 | as.factor(treatment) + as.factor(placebo) + turnout2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,])
summary(cace_turnout)
#Column 4
cace_turnout1 <- ivreg(turnout2017 ~ treat_rate + turnout2015 | as.factor(treatment) + turnout2015, data = FieldExperimentData[!FieldExperimentData$letter_condition == 0,])
summary(cace_turnout1)

######################
# Online Appendix L  #
######################

summary(neighbor_model <- lm(proindependence2017 ~ neighbor_treated + proindependence2015, data = FieldExperimentData[FieldExperimentData$letter_condition != 0 | 
                                                                                                                    FieldExperimentData$neighbor_treated != 0,]))




